####################################################################
# Figure 7
#####################################################################

#####################################################################
# Preliminaries
#####################################################################

rm(list = ls())

library(scinference)
library(limSolve)

#####################################################################
# Functions
#####################################################################

source("code/auxiliary scripts/common_functions.R")

#####################################################################
# Calibration DGPs
#####################################################################

source("code/auxiliary scripts/calibration_dgps.R")

#####################################################################
# Simulations
#####################################################################

set.seed(12345)

nreps <- 10000

tau_grid_g <- c(seq(0,0.4,0.05))

all_DGPs_g <- array(NA,dim = c(length(tau_grid_g),4,9))

for (DGP in 1:9){

  results_overall_g <- matrix(NA,length(tau_grid_g),4)

  for (t in 1:length(tau_grid_g)){
    
    results_g <- matrix(NA,nreps,4)
    
    for (r in 1:nreps){
      results_temp_g1 <- sim_one_sample(DGP,T0=30,T1=T1_app,J=J_app,K=4,Lambda,rho_u,var_u,var_factors,rho_vec,var_epsl_vec,w0_sc,tau_alternative=tau_grid_g[t],li=FALSE,sdid=FALSE,permtest=TRUE)
      results_temp_g2 <- sim_one_sample(DGP,T0=150,T1=T1_app,J=J_app,K=8,Lambda,rho_u,var_u,var_factors,rho_vec,var_epsl_vec,w0_sc,tau_alternative=tau_grid_g[t],li=FALSE,sdid=FALSE,permtest=TRUE)
      results_g[r,1:2]  <-  1-results_temp_g1$cov[c(1,4)]
      results_g[r,3:4]  <-  1-results_temp_g2$cov[c(1,4)]
    }
    
    results_overall_g[t,] <- colMeans(results_g)
  
  }

  all_DGPs_g[,,DGP] <- results_overall_g

  # print(DGP)
  
}

#####################################################################
# Power curves
#####################################################################

for (DGP in 1:9){
  
  results <- all_DGPs_g[,,DGP]

  graphics.off()
  pdf(paste("graphics/power_comparison_",DGP,".pdf",sep=""),pointsize=16,width=8.0,height=6.0)
  plot(NULL, xlab="Alternative", ylab="Power", main=paste("DGP",DGP), col="blue", pch=".", xlim = c(0,max(tau_grid_g)), ylim=c(0,1))
  lines(tau_grid_g,results[,1],col="black",lwd=4, lty=1,type="b",pch=1)
  lines(tau_grid_g,results[,2],col="gray",lwd=4, lty=1,type="b",pch=3)
  lines(tau_grid_g,results[,3],col="black",lwd=4, lty=2,type="b",pch=1)
  lines(tau_grid_g,results[,4],col="gray",lwd=4, lty=2,type="b",pch=3)
  abline(h=0.1,lty=1,lwd=1,col="gray")
  if (DGP==1){ legend("right",legend=c("t-Test, small","Permutation test, small","t-Test, large","Permutation test, large"), seg.len=2, col=c("black","gray","black","gray"),fill=NA,border=NA, lty=c(1,1,2,2),lwd=c(2,2,2,2), pch=c(1,3,1,3), merge=T,bty="n")}
  dev.off()

}
